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A METHOD FOR DYNAMIC CHARACTERIZATION OF 
DENSITY FIELDS IN A COMPOUND STRUCTURE 

5 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

The present invention relates to the field of characterizing the properties of a 
1 0 field, particularly to the characterization of properties within a field relating to mass 
transfer or material movement. The characterization of field properties is related to 
fields such as localized absorption/adsorption of materials to surfaces or volumes, 
especially in the delivery of bioactive materials to in vivo tissue during therapeutic 
treatments. 

15 

2. Background of the Art 

It is a routine procedure in image processing to "segment" a two-dimensional 
(2D) or three-dimensional (3D) image. This can be performed, for example, on 1) a 
photograph of a slice of biological material, or 2) a scan by computer-aided 

20 tomography (CAT), magnetic resonance imaging (MRI), sonogram, or seismology 
recordation systems, where the data resembles a stack of slices. Many ways of 
gathering such data exist in various areas of technology, such as, for example, 
scanning data from individual components in an image to provide distinct layers 
with related data (e.g., only single colors, only ranges of optical density, etc.). 

25 The image provided in such segmenting comprises a rectangular array of 

points with specified intensity, density or color values, the points being called 
"pixels" in 2D and "voxels" in 3D. The objective of the segmentation procedure is 
to classify these points into classes of materials so that information can be read into 
or from the image. In a specialized technological field, such as in medical 

30 applications, for example, "bone", "muscle", "blood", or "tumor" may be assigned to 
different points and ultimately areas of the image. In geology, the classes might be 
"granite", "shale", "water", etc., and so on for other sciences. The typical starting 
point for classification of the image is the value attached to a specific point. For 
example, bone absorbs X-rays more than any other tissue, so that one can fix a 



threshold value T according to the scaling of a particular Computed Axial 
Tomography (CAT) image and use "opacity greater than T" as a test for attaching 
the label "bone" to a point. A "segment" is then either a) the class of points so 
classified, or b) a connected component of this class. For example, it can be useful 
5 in medical applications to recognize the left kidney and the right kidney as separate 
segments of the body. Various methods also exist for correcting errors and noise in 
the image. For example, an isolated voxel with a value associated with muscle, 
surrounded by voxels classifiable as bone, should usually be labeled as bone also. 

There exists also a large range of software to assist humans in segmentation, 

10 allowing them to draw curves around a region so that the regions can be recognized 
as separate. For example, where a region presses closely against another region 
with similarly associated values (such as a bone touching a bone), single values 
alone would be inadequate for classifying points. Despite such software assistance, 
this region recognition technology is labor intensive, especially in 3D where the 

15 classification is usually done slice by slice. Therefore, improved automation is a 
continuing area of vigorous research. We assume herein that it is desirable that 
segmentation includes not only a labeling of grid points by segments, but also that 
some descriptive geometric information such as bounding areas such as its bounding 
box (the smallest rectangular region with sides parallel to the x, y and z axes that 

20 contains every point of it) or a list of boxes in some scheme such as an octree 
division of the grid is provided, so that between all of these components of the 
segmentation labeling system, all of the points of the segment are included and 
identified (classified). Such additional data allow iteration over a shorter sequence 
of grid points than the whole set, when points outside the segment will contribute 

25 nothing to the intended result of the classification within the segmentation. Where 
such data are not delivered with the segmentation, they can be rapidly constructed 
for the entire segment list by a single run-through of the grid. Therefore, we 
assume that they are available. 

Once 3D segmentation is achieved, it allows volwnetry, the defining and/or 

30 measuring of the volume of a segment, which is proportional to the number of 

voxels labeled as belonging to it. This can be important either in static images, such 
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as a geological or body scan where it answers questions like "How much ore is 
present?" or "Is the liver enlarged?" or in a sequence of 3D images, such as a record 
of a beating heart, where the change in volume of the left ventricle can assist in 
quantifying what fraction of the blood is expelled into the aorta, and thus how 
5 effectively the heart is functioning as a pump. As a point is either In or Out of a 
segment (it is or is not blood, muscle, granite, etc.), volumetry techniques have 
developed for this two-valued context. The present invention addresses the 
computation of quantities appropriate to, and by methods appropriate to, 
continuous- valued variables such as a density or concentration. 



A method that supplies a unified suite of quantification functionality for 
density functions defined on or defining a three-dimensional space, which may 
optionally vary in time, includes at least two services assisting in the definition. 
1 5 These services may include in a dynamic modality: 
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SUMMARY OF THE INVENTION 



(a) Computation of the volume of the region where the density lies 
above or below a specified threshold, or between two specified 



values. 



20 



(b) Computation of the integral of the density (that is, determining a total 
amount of material within the region). 



(c) Estimation of the rate of change of the density with respect to time, 
optionally restricted to a specified region. 



25 



(d) Estimation of the local or global failure of conservation represented 
by changes with time in the density, whether or not the method is 
given an implemented transport model. 



(e) Estimation of the local or global rate at which material with a 
changing density is passing through a specified surface (e.g., 
boundary), whether or not the method is provided with an 
implemented transport model. 



30 



(f) Separation of the density of a material with a changing density, given 
an implemented transport model, into "free" and "bound" densities. 
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The present invention does not address the static segmentation of an image 
into classes, as described above. Rather, the invention assumes such a classification 
to have been completed, if needed. The present method addresses the dynamic 
5 volumetry of a first material diflusing or being otherwise transported through 
another substance (e.g., the medium or environment through which the diffusing 
material is passing), that substance through which materials are being transported 
may be, by way of non-limiting examples only, a compound structure such as a 
brain, terrain, waterway, habitat or mountain. We refer to the latter substance, 

10 whether simple or compound, as the substrate. Various methods have recently 

made it possible for scanning systems to estimate the densities c of such a substance 
at a grid of positions (x, y, z), such as labeling drugs for MRI visibility (e.g., U.S. 
Patent Nos. 6,061,587 and 6,026,316) or attaching green fluorescence protein (GFP) 
that makes the material glow under suitable conditions. Equally, such an array c (x, 

15 y, z) of scalar values c (x, y, z) may arise from a simulation of a transport process, 
often for purposes of simulation, as in [Copending U.S. patent Application bearing 
attorney's docket number 723.032US1, titled "A METHOD AND APPARATUS 
FOR TARGETING MATERIAL DELIVERY TO TISSUE" baring U.S. Serial 
No. 09/566,478, filed May 8, 2000]. It may also come from the evaluation at grid 

20 points of an algorithmically specified function: for example, the most concise 

specification of "the ball of radius R centered on (X Z Z)" is as the set of points (x, 
y z) where c (x, y z) = {{x - X) 2 + (y-Y) 2 + (z- Z) 2 ) /R 2 does not exceed L This is 
an implicit specification, defined by a testable criterion rather than by a finite or 
infinite list of points. With more complicated functions, c, implicit functions are 

25 widely used in computer-aided design (CAD). The present invention may be 

applied to such c also, allowing (for instance) efficient calculation in CAD of both 
volumes and - where constant or variable mass density is known - weight. 
Throughout this document, the symbol c refers to a scalar quantity (number) for 
which possibly different values exist at a plurality points in space and where 

30 appropriate at a plurality of times. (That is, c is a function of position or of position 
and time.) Sometimes the plurality of points at which values are available is 
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continuous, sometimes it is restricted to a grid of points indexed by integer triples 
(ij\k) with indices lying in a specified range, at a fixed time or at a finite plurality 
of times. (This may arise by storage of an array of c values at such a grid found by 
physical scanning, or by evaluation at grid points of a function defined on arbitrary 
5 points, at the said fixed time or plurality of times.) Where it is necessary to 
emphasize that c is available only at grid points, with values stored in an array 
(which then formally constitutes the function produced by restricting c to the set of 
those grid points), we distinguish it as c g . If c is a temperature in tissue subjected to 
electromagnetic fields in an imaging system, the present invention applies to such 

10 questions as 6t what fraction of the hippocampus is heated to a temperature above 
38°C?" If c (x, y, z) is the heating rate at (x, y, z), we may wish to compute total 
heat input, and so on. 

For generality, therefore, we refer to c (x, y, z) or c (x, y, z, f) as a density, 
and to the material, heat, rate or other real or formal quantity that c describes as a 

15 superstrate. 

The present invention applies to all situations in which such a grid c g of c (x, 
y, z) information can be estimated during dynamic variations in the superstrate. The 
interest is particularly in the case in which the superstrate is moving, so that 
repeated scans provide a sequence of changing value grids, that is a dynamic image 

20 of the material (condition, property, etc.) and its relationship within the superstrate. 
Note that the signal strength found by the material detection system is not always 
linearly proportional to the density, but as long as the measured quantity changes 
monotonically with the density and has a consistent value at each density, c can be 
deduced from it. We assume such pre-processing to have been done, and work with 

25 a grid c g of density estimates c. In certain cases the scanning system can distinguish 
between bound states of the superstrate, where the superstrate is attached to a 
particular point of the substrate, and free states, where the superstrate can move 
through the substrate. In such cases, we distinguish the bound density c b from the 
free density cf. Often a superstrate has a desired therapeutic effect, or an undesired 

30 effect, only when its density lies between certain limits. It is thus useful to quantify 
the volume of substrate in which c or c b or <f is within such a range. 
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Given an estimate c g of the density c, or c b and </, the present invention 
computes (/) the overall volume of the region in which c or c b or <f lies between 
user-specified limits Q ow and Chigh, (w) the total amount of superstate, (Hi) the rate 
at each point at which the superstate is being created ex nihilo or from other 

5 distributed entities not contributing to the estimates c, c b , or <f 9 less the rate at which 
it is being destroyed, with or without production of other entities, (iv) the total rate 
at which the superstate is being created, less the rate at which it is being destroyed, 
(v) where the scanning technology permits, the separate totals of free and bound 
superstrate, (vi) the rate at each point at which the superstate is being bound, less 

1 0 the rate at which it is being freed at that point, (vii) the overall rate of binding less 
freeing of the superstrate, (viz/) the rate at which the superstate is leaving or 
entering the scanned region through each point on its surface, and/or (ix) the overall 
rate at which the superstrate is leaving or entering the scanned region. The present 
method further computes for any user-selected region, S 9 the corresponding 

1 5 restricted quantities (I) the overall volume of the region within S in which the 
density lies between user-specified limits Ci ow and Chigh, (II) the total amount of 
superstrate with S 9 (III) where the scanning technology permits, the separate totals 
of free and bound superstrate within S 9 (IV) the total rate at which within S the 
superstrate is being created ex nihilo or from other materials not contributing to the 

20 estimate c, less the rate at which it is being destroyed, (V) the overall rate within S 
of binding less freeing of the superstrate, (VI) the rate at which the superstrate is 
leaving or entering S through each point on the surface of S 9 and/or (VII) the overall 
rate at which the superstrate is leaving or entering S. 

Given a segmentation of the scan region that is spatially registered with the 

25 density scan (and in the case of a sequence of segmented scans of a moving 

structure or superstrate or sub-component therein, temporally registered also), the 
region S can be a segment such as a tissue T 9 allowing the system to report such 
quantities as "total material in the hippocampus". 

30 BRIEF DESCRIPTION OF THE DRAWINGS 

Drawing 1 shows a sum of values serving as an approximate integral. 
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Drawing 2 illustrates that a piecewise linear fit to data points gives a 
different estimated integral than a spline fit. 

Drawing 3 illustrates the error in estimation of content (area or volume) by 
counting grid points. 

5 Drawing 4 illustrates a method using regions tipped by pieces of boundary. 

Drawing 5 shows the distinct ways (up to rotation and in/out exchange) that 
a volume boundary can meet the edges of a volume element. 

Drawing 6 shows the data and function members of the object class that 
delivers the time-independent functionality of the invention. 
10 Drawing 7 shows the data and function members of the object class that 

delivers the time-dependent functionality of the invention. 

Drawing 8 shows the logical flow of the computation of the volume enclosed 
by a set of triangles. 

Drawing 9 shows the logical flow of a coarse, fast estimate of the volume in 
1 5 which the density lies between limits C\ and Cj. 

Drawing 10 shows the logical flow of a more precise estimate of the volume 
in which the density lies above or below a threshold C 

Drawing 1 1 shows the logical flow of a version of the corresponding 
estimate of the volume within a segment S in which the density lies above or below 
20 a threshold C 

Drawing 12 shows the use of the logic in Drawing 10 or 1 1 to give the 
corresponding estimate of the volume within a segment S in which the density lies 
between limits C\ and C2. 

Drawing 13 shows the logical flow of the computation of a free density from 
25 a free-plus-bound density and a transport model. 

Drawing 14 shows the logical flow of the calculation of net efflux through a 
triangulated surface. 

DETAILED DESCRIPTION OF THE INVENTION 

30 A method supplies a dynamic vector map of properties within a region, or as 

otherwise stated, a unified suite of quantification functionality for property 
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functions, such as density functions, conduction functions (e.g., thermal conduction, 
electrical conduction, atomic or subatomic mass conduction, macromolecular mass 
conduction), defined on, defined in or defining a three-dimensional space, which 
functions may optionally vary in time. The method should include at least two 
5 services assisting in the definition of the map or suite. These services may include 
as non-limiting examples, in a dynamic modality: 

(a) Computation of the volume of the region where the density lies 
above or below a specified threshold, or between two specified 
values. This type of service assists in providing a basis for 

10 characterizing a point or region where a property (e.g., density) can 

be a part of the definition of the substrate. 

(b) Computation of the integral of the density (that is, determining a total 
amount of material within the region). This provides a basis for 
establishing a foundation for comparison of changes within the 

15 region, as a starting point for determining changes must begin with 

an initial or ground state value. 

(c) Estimation of the rate of change of the density with respect to time, 
optionally restricted to a specified region. Rate functions are critical 
in determining dynamic properties for mapping. The rate may be 

20 with respect to concentration changes, conductivity rates, 

temperature changes, optical density changes, viscosity changes, or 
any other observable property. 

(d) Estimation of the local or global failure of conservation represented 
by changes with time in the density, whether or not the method is 

25 given an implemented transport model. This service can identify 

changes indicating events wherein the total superstrate amount does 
not remain constant within the region or volume. This could be 
evidence of mass transfer out of the region, destruction or creation of 
material within the region (e.g., expression or metabolism), and other 

30 events that alter the integrated totality of measured substrate within 

the region. 
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(e) Estimation of the local or global rate at which material with a 
changing density is passing through a specified surface (e.g., 
boundary), whether or not the method is provided with an 
implemented transport model. This service is an indication of 

5 transfer rates between segments within the region, as between solid 

surfaces and moving liquid mass in contact with the surface. 

(f) Separation of the density of a material with a changing density, given 
an implemented transport model, into "free" and "bound" densities. 



10 We assume that a scanning system (hardware and software), simulation, or 

function sampling delivers a grid c g of density estimates c$k at points rectangularly 
arranged in space, for integer values 0 </ < /, 0 <j < J and 0 < k < K, with positions 
(x, y, z) = (ih X) jh y , kh z ) in a suitably chosen system of rectangular coordinates. The 
ranges of the integers i, j and k may begin other than with 0, and one skilled in the 

15 art may easily extend the application of the present invention to the case where the 
region sampled is not rectangular, but since a rectangular array of points in a 
cuboidal region is overwhelmingly typical for the presentation of scan data, we use 
this case for exposition. Similarly, the points may be arrayed in the face-centered 
cubic lattice arrangement typical of the centers of stacked cannonballs, or according 

20 to any other regular scheme, without materially altering the application of the 
present invention. Other equivalent volumetric formats for analysis may also be 
used, such as spherical arrays or regions, pyramidal arrays or regions, or any other 
geometric formats or mathematical formats. There does not at this time appear to be 
any reason for assuming significant benefits in the practice of the invention for 

25 using formats other than cubic, rectangular or spherical regions, but the others may 
be used as a form of equivalent which might prove to be desirable because of the 
shape or certain three-dimensional volumes that are being considered in specific 
applications of the process to a patient or treatment. Mixtures of volumetric formats 
and regions may, of course, also be performed and used in the practice of the 

30 invention, especially where tighter fit (higher accuracy) is desired and the software 
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can select from amongst different shapes and formats to determine best 
approximations of volumetric fits. 



Superstrate quantification 

We denote the region sampled as R, and assume that the amount (volume, 
mass, conductivity, adsorptivity, energy, or any other property or characteristic) of 
superstrate within R is the integral over R of the density c, of which the data grid c g 
provides estimates Cyk only at grid points. A simple estimate of this integral is given 
by 



10 MA2>* (!) 

i=0 
J=0 
k=0 

since h x h y h z is the volume in (x, y, z)-space of the cuboidal volume element of points 
nearer to a given (x, y, z) = ih x , jh y , kh z ) than to its neighbors in the grid. In one 
dimension, this corresponds to estimating the integral of a function f(x) by a sum of 
the values c t =f(ih), 

15 (2) 

which assumes that/ has a bar-chart graph as shown in Drawing 1, and sums the 
areas of the bars. An alternative approximation / to / is to interpolate linearly 
between the sample points and find the area 



h 7-1 

— 2( c * + c i ) where / = / + 1 (3) 

2 i=o 

20 under the resulting polygonal curve, but the result is identical except that (2) 

extrapolates for half a step beyond the end points. If the end contributions are made 
proportional to the length element fraction that is not outside the end points, the 
result coincides with (3). Similarly, if we replace (1) by trilinear interpolation 
within each rectangular volume element whose edges are axis-direction steps in the 

25 sample array, the resulting estimate: 
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£-1 

KKK S 

— — X {c lfi + c- ljk + c h]k + c - + + c fk + c- ifi + c m ) (4) 

O i=0 

agrees with (1) if the latter is adjusted to reflect the half, quarter and eighth volumes 
(within the grid) of the elements meeting the boundary. 

However, if an f(x) known only by its values at step points is approximated 

5 by a smoother curve / such as a cubic spline (Drawing 2), the area under / (and 
hence the estimated total for J) is slightly different. This result is similar for 
smoother fits in three dimensions. This is pointed out only to observe that there are 
multiple means of estimating the total of a sampled density, of which the present 
invention may use any method for sampling density of the superstrate or the 
10 substrate. For simplicity, a preferred implementation is the use of (1), whose 

difference from other estimates is typically less than the uncertainty resulting from 
noise in the sampling process. The present invention is easily applied by one skilled 
in the art to the use of other estimation methods. 



15 Volume quantification 

Defining a range function p as 

^\\#p<x<q 
p,q | 0 otherwise 

^ f lif jc<? 

^ >q [0 otherwise 

or/? (*) = { lifp ' x (7 ) 
p * [0 otherwise ' 

20 this provides a "density within range" volume 

P( c ( x >y> z))^xdydz (8) 

which may be estimated in various ways for the purposes of the present invention. 
The simplest estimate is performed by counting the array points (x, y, z) = (ih x , jh y , 
khz) for which r(c$) - 1, equivalently summing the function p (c). However, as 
25 Drawing 3 illustrates for the two-dimensional case of computing area rather than 



11 



volume, this introduces unnecessary inaccuracy. The area ar 2 ~ 3.1416 of the set 
{(*> 7) | x 2 + / < 1 } is poorly estimated as 3.17 by counting the 3 17 squares whose 
centers lie inside the circle. Several more accurate methods may be used. For some 
specific fields of use, this degree of inaccuracy would be acceptable, while with 
others it might not be tolerable. 

In one alternative method, one first creates a piecewise linear description of 
the boundary of the range: in two dimensions a polygonal curve, in three dimensions 
a polyhedral surface. Drawing 4 shows the circular edge of the set {(x, y) \ x 2 + 
y*< 1 } approximated by an oriented loop 401 of twenty-four straight edges. 
Orienting this loop anticlockwise, the areas of strips joining each strip to a side 410 
of the region R sum to the area inside the loop 401 , if opposite signs are applied to 
strips like 420, where the edge 421 projects to an upward step 422 on the side 410, 
and strips like 430, where the edge 431 projects to a downward step 432. 
Analogously, suppose that the surface {(x, y, z) \f(x, y, z) = C} is approximated by a 
closed set S of triangles, as may be done by the Marching Cubes algorithm (in 
Lorensen, W. E. and Cline, H. E., SYSTEM AND METHOD FOR THE DISPLAY OF 

SURFACE STRUCTURES CONTAINED WITHIN THE INTERIOR REGION OF A SOLID BODY, 

US Patent 4,710,876, 1987), the Skeleton Climbing algorithm (T. Poston, T-T. 
Wong, and P-A. Heng, "Multidirection Isosurface Extraction with Adaptive 
Skeleton Climbing, Computer Graphics Forum, 17, 3, September 1998, pp. 137- 148; 
and T. Poston, H.T. Nguyen, P-A. Heng, and T-T Wong, "Skeleton Climbing: Fast 
Isosurfaces with fewer Triangles," Proceedings of Pacific Graphics '97, Seoul 
Korea, October 1997, pp. 117-126., or otherwise. If each triangular face a e Z has 
vertices 

v 0 O) = (* 0 OX ^0 OX z o O)) 

v 1 (o-) = (x 1 ( < t),j; 1 (o-),z 1 ( ( t)) (9) 
v 2 (<t) = (x 2 (a), y 2 (a), z 2 (a)) 

so ordered that they appear anticlockwise when viewed from outside the surface, we 
may write the volume inside the surface as a signed sum of volumes of slant-topped 
prisms 
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£ (( Xi (a) - x 0 {(j)) {y 2 (a) - y 0 (<j)) - (x 2 (o) - x Q (<jj){y x (a) - y 0 (<7))){z 0 (<x) + z x (d) + z 2 (a)) 

aeZ ______ 

6 

(10) 

Note that if the set S meets the higher-z face of R, the area of the part of R within the 
volume (according to the approximation I of its boundary) must be included for a 
5 correct result. Topologically triangles could be included within any face of R met 
by the in-range volume, for a closed surface, but triangle within the sides and 
bottom contribute nothing to the volume sum (10) and may thus be omitted from the 
calculation. 

The method above is an efficient volume computation where such a set 2 has 
10 already been constructed (for example, in order graphically to display the region 
where c is within the range of interest), but constructing £ solely for this purpose 
represents unnecessary computation. In a further alternative method, we classify 
volume elements with corners on the sample grid according to the way their vertices 
meet the in-range and out-of-range values of c, and estimate the volumes of the 
1 5 within-range sets within them. 

Drawing 5 shows different ways that the vertices of a volume element may 
be within range (shown as a full circle) or out of range (shown as an empty circle), 
with a corresponding simple polygonal curve in which the surface of the in-range 
set can meet the faces of the element, consistent with the data. This is similar to the 
20 classification in Lorensen, W. E. and Cline, H. E., "SYSTEM AND METHOD FOR THE 

DISPLAY OF SURFACE STRUCTURES CONTAINED WITHIN THE INTERIOR REGION OF A 

SOLID BODY," US Patent 4,710,876, 1987, hereinafter referred to as the Marching 
Cubes algorithm, but differs where a face is ambiguous, with exactly two corners 
within range and these are diagonally opposed. In an ambiguous face, we so choose 

25 the curve as to separate the in-range and not the out-of-range corners. Other 
consistent rules can be followed within the logic of the present invention. (For 
example, the opposite rule, or the rule that the curve length should be as short as 
possible may be used.) Every one of the 256 ways that the eight corners may 
between them be in or out of range is represented in Drawing 5, by a suitable 

30 rotation and/or a reversal of in-range and out-of-range, with reversal not permitted 
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where there is an ambiguous face. (Permitting reversal in such cases would lead the 
original version of the Marching Cubes algorithm, as first described by Lorensen 
and Cline in the above-cited patent, to create triangle sets that are not as there 
claimed, separating surfaces between the in-range and out-of-range sample points.) 
5 It is not efficient in an implementation to determine to which of these cases 

an element may be reduced, and to rotate and/or reverse it to standard form, but it is 
more convenient here to show 17 cases rather than 256* Our preferred 
implementation uses a lookup table to find an appropriate transformed and 
implemented version of one of the formulas below, according to which of the 256 

1 0 cases it represents. 

The vertices of the polygonal curves shown in Drawing 5 lie on the element 
edges that cross between in-range and out-of-range sample points. Their position 
may most conveniently be estimated by linear interpolation between the sample 
values at the edge ends, but more complex procedures (such as those suggested in 

15 various papers elaborating marching Cubes) are usable within the framework of the 
present invention. For each such vertex, Drawing 5 labels the distance from one or 
the other neighboring sample point on the volume element edge. 

There is no unique choice of surface spanning the polygonal curves in 
Drawing 5, but simple choices constructed from flat triangles and quadric 

20 quadrilaterals lead to the following formulas for the volumes contained by the 
surface within the volume element, all to be multiplied by the factor h x hyh z l\2. 
Where a case not shown can rotate to a case in Drawing 5, the same formula 
applies: where reversal of in-range and out-of-range points is needed, the formula 
below should be subtracted from the number 12. 
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0 

labc 

(2(bc + de) + be + dc) 
2(abc + def) 
2(ab-abf + b + e- fb~ fe + bc + ce + 2dec-2de) 
2(abc + def) 
2(d-ae($ + b) + \0(c + e)-9ac) 
2(bc + de + fgh) + be + dc 
2(6-hg(8 + c)-b-d-c-bc + bg(c + f -l)-e- / + 4g + 3h + f (d - g) + gc + 2hc + he) 

2(abc + def + ghi) 
2(5-b + ba(g-l)-g-i + ib-hg-ih + ij(h-l)-def) 
3(a + d + b + c) 
2(5-b + ba{g-\)-g-i + ib-hg-ih + ij(h-\)) 
2(bc + de + fg + ih) + be + dc + fh + gi 
2(ac + af + / + e-ef +b-3be-ba + a + de) 
2(d-ae$ + b) + \Q(c + e)-9ac + ghi) 
2 (abc + def + ghi + jkl) 

(ID 

Where the desired volume is that of the points where c is between two limits C\ < 
C2, this algorithm may be used to find the volume below C\ and subtract it from the 
5 volume below C2, either volume element by volume element or as a total, depending 
on details of architecture and efficiency. 

Transition Quantification 

Superstrate density at a point (x, y, z) may change because superstate 

10 moves, or because some transformation creates or annihilates the superstrate as a 
detectable entity with the scanning method currently in use. If the superstrate is a 
material inserted into the body with a "contrast agent" label whose detectability is 
constant (such as an atom that heavily absorbs X-rays, or a gadolinium atom visible 
to appropriately tuned MR imaging), only movement will be detectable. Certain 

1 5 labels will change frequency when their chemical configuration or chemical 
structure alters (as with F-19, Fluorine- 19), and changes in MR response 
frequencies could then be used to assist in determination of changes in the state or 
volume of a specific compound, based upon the changing responses due to changes 
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in chemical state of compounds with the label. (See, for example, the techniques 
used in copending ILS. Patent Application, titled "Imaging Methods For Visualizing 
Implkanted Living Cells," Attorney Docket No. MAL 500.003US1, Inventors 
Michael Moseley et al., filed on June 28, 2000, bearing U.S. Serial No. 09/606,137, 
5 which imaging techniques and methods, particularly by Magnetic Resonance 

Imaging, are incorporated herein by reference.) For example, many applications of 
superstate volumetry involve substances naturally created or destroyed by 
biological processes, and detectable by intrinsic properties. For example, the 
magnetic resonance response of an atom depends on its chemical situation. A 

1 0 hydrogen nucleus fixed in a large molecule is less free to rotate than one in a free 
water molecule, and this difference is one of the main contrast sources in MR 
imaging. Bound water, as in certain crystals, is distinguishable from free water. Ice 
is similarly distinguishable from free water, so that MR offers the possibility of 
sensors appropriate to control of cryogenics, managing the freezing of biological 

15 tissue to minimize damage. 

With appropriately tuned MR pulses, such detectability may be made highly 
specific for nuclei in a particular molecular setting, thus making the MR scan a 
direct detector for a particular molecule in a particular state. 

Given only a grid c g at discrete points (ij\k) of samples or local estimates of 

20 a changing record c (x, y, z, t) at a sequence of times t 9 it is impossible to distinguish 
the net rate C (x, y z, t) of creation (creation less destruction) of the superstate with 
density c from the effects of superstate transport. (In particular, it is logically 
allowable that all increments and decrements of c are due to local creation and 
destruction.) However, if the transport properties of the superstate are known, we 

25 may write 
dc 

—(x 9 y 9 z 9 t) = (netarrival(x,.y,z,0)- C{x 9 y 9 z 9 i). (12) 
at 

In the case of pure diffusion, for example, the arrival term is D (jt, y, z) V 2 c(x, y, z, 
t) 9 where the field D gives the diffusion tensor at each point. (In simple cases, this 
may reduce to a single constant, independent of position (x, y> z) and the direction of 
30 diffusion.) A copending patent application (see copending U.S. Patent Application 
Serial No. 09/566,478, bearing attorneys docket no. SLWK 723.032US1, titled "A 



METHOD AND APPARATUS FOR TARGETING MATERIAL DELIVERY 
TO TISSUE" filed on May 8, 2000) describes how such terms may be established 
for a given substrate, with a model-dependent distinction between transport and 
destruction. (For example, when a drug diffusing in brain tissue may pass into veins 

5 and is removed by the blood flow, this may be treated as a "sink" rather than as 
explicitly modeled transport.) Where the arrival term is thus computable, the left 
hand side of (12) may be estimated from data, in general using some smoothing 
with respect to space and time since differences like c^ijjc /next)- c g (y,£, W) are 
more strongly affected by noise than are the estimates of c itself. The net creation 

1 0 rate C thus becomes a known scalar field such as c. The procedures described 
above can be applied to this, finding the total rate of creation by estimating its 
integral, the volume of the region in which creation occurs at a rate above some 
threshold, and so on. (In this case our "density" is a rate, and our "superstate" is a 
process.) 

15 In some applications, it is important to distinguish the bound density c b from 

the free density <f. In some cases, the imaging technique may be tuned to this 
difference (for example, the binding of a molecule may change its nuclear resonance 
properties), so that c b and <f are directly detectable, and their volumetry may be 
handled separately by the methods above. In cases where superstrate is conserved, 
20 there is an alternative means of separation. 

This alternative can be performed as for the transport dynamics applies only 
to the free density, so that 
dc f 

— - (x 9 y 9 z 9 t) = (net arrival (x 9 y 9 z, f )) - binding rate (13) 
at 

(x 9 y 9 z 9 t) = binding rate (14) 



d^_ 
dt 



dc d(c f + c b ) 

25 —{x 9 y 9 z 9 t) = (x 9 y 9 z 9 t) = (netarrival(x,jF,z,0) (15) 

at at 

Taking the case of diffusion transport as an example, the net arrival rate is D (x,y,z) 

VV (x,y,z,t) 9 so that for the transport operator T= D (x,y,z)V 2 : 

f dc 

Tc J (x 9 y 9 z 9 t) = —(x 9 y 9 z 9 t) (16) 
at 
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where the right hand side is known by direct or indirect measurement. Given 
boundary conditions (such as a boundary on which c is zero, to which the 
superstate has not yet diffused), knowledge of V V implies a unique solution for J . 
The same is true for more general transport operators T, such as many of those 

5 arising in (the above mentioned attorney 5 s docket number SLWK 723 .032US 1 , 
titled "A METHOD AND APPARATUS FOR TARGETING MATERIAL 
DELIVERY TO TISSUE" U.S. Serial No. 09/566,478, filed May 8, 2000). The 
data and logical structure information necessary to invert Tare aspects of the 
transport model For the purposes of the present invention, it is thus assumed that 

1 0 the transport model includes a function which, given access to a dc/dt value array on 
whose boundary c is negligible, fills another array with the solution of (16) for cf. 

Thus the values are available for c b =(f-c, and changes may be tracked in c b 
to give the binding rate. All these quantities may be separately displayed and 
integrated by the methods above. For computation, we replace the various c, </, c h 

15 by their corresponding grid estimates, and derivatives by the corresponding finite 
differences. 

Segmented Superstate Volumetrv 

Given a segmentation of the substrate (for instance into bone, white matter, 

20 gray matter, etc.), any of the above computations may be restricted to a user- 
selected segment S. In such a case the algorithm applies unaltered to each volume 
element that has all eight corners in S, and give 0 where all are outside S. In a 
doubly-boundary volume element E that contains a non-zero number of corners 
inside and corners outside S 9 and of corners above and below threshold, it is 

25 necessary to estimate the volume within EnS that is on the desired side of the 

threshold. This may be done by various methods, according to the precision desired 
(since doubly-boundary elements are relatively few, the percentage impact on the 
total volume estimate is proportionally small). Most coarsely, one may multiply the 
estimate above by n/S, where n is the number of corners of E that are in S. More 

30 accurately, one may estimate the fraction of surface area of E that is in S, and use 
this as a multiplier. More accurately again, one may use 2D geometry in each face 



18 



to estimate the fraction of surface area of E that is both in E and to the desired side 
of the threshold, and return h x h y h z times this fraction. (This is equivalent to a 
surface model for each surface area that places a vertex in the middle of E, and 
creates a surface triangle for each straight edge in a face.) Most accurately, one may 

5 construct a triangulated surface independently for each surface area, add vertices 
and edges where the edges and faces of these surfaces intersect, add triangles to fill 
in the resulting network of edges, and compute the volume enclosed. This makes 
the code more complex, but for a small number of doubly-boundary elements it does 
not greatly increase computation time except in an implementation where all 

10 elements treated in parallel in a given step must be completed before the next 
parallel step begins, where processors which have quickly resolved simple cases 
would have to wait for the completion of doubly-boundary element computations. 

Where the substrate itself is in motion (for example, a heart or lung), the 
segment S need not be a static region in (x,#z)-space, but may be a region S (t) that 

1 5 depends on time. The implementation thus requires a representation of moving 
regions. 

Rates of change 

Where the superstrate density is a function c (x,y,z,t) of both position and 
20 time, it is often important to report total changes and their rates. In a setting with 
neither noise nor quantization error, the integral over a region S of the rate of change 
of density is identical to the rate of change of the integral of density over S. In 
practice, since rate of change estimation is noise sensitive, this identity depends on 
the estimation and smoothing methods used. In the present invention, for reasons of 
25 numerical efficiency, the preferred implementation is first to compute estimates of 
"density within range" volume and of superstrate total quantities at the available 
times /, and then to estimate their rates of change by fitting smooth functions of t to 
the resulting sequences of scalar values. For numerical integration, these functions 
are typically evaluated at a mesh of grid points (ij,k). 

30 
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Software structure 

In one embodiment, the present invention does not constitute stand-alone 
software with an end-user interface. It can provide important functionality to a 
volume-exploration program, which also provides rendering routines for display of 

5 surfaces and volumes, and a set of mousepad or reach-in widgets for the user to 
control what is shown. It provides essential functionality to control software for 
image-guided delivery of drugs or living cells, as in U.S* Patent Nos. 6,061,587 and 
6,026,316, with needs in its display functions and user interface that overlap with, 
but are distinct from those of a general volume-exploration program. It could 

10 provide essential functionality to (as a hypothetical example) software able to 
address time-sequenced three-dimensional scans of a hundred rats breathing 
contrast-labelled tobacco smoke — thousands of volume scans, in all — and return 
statistics on the depth of penetration at particular smoke particle densities, with no 
user viewing of individual images. 

1 5 In another embodiment, the method is implemented as a software library, 

with an application-programming interface (API) that may be used by any 
application, such as those just mentioned. The fundamental software object is the 
"Quantifier" class 600 shown in Drawing 6, which has data members 610 and 
function members 630. Any valid Quantifier 600 must have specified dimensions 

20 61 1 J, K for the array which contains the superstate density values, and know the 
type 612 of number (integer, floating point, etc.) which is stored there, and the 
address 613 of the array in which superstate densities are stored. These members 
may be defined for it by member-setting functions 640, which may take a variety of 
arguments, acted on in different ways. (For example, if the function that sets the 

25 superstate array receives a file name or Uniform Resource Name as argument, it 
will allocate memory and read the data into it. If it receives the address of data 
already in RAM and accessible for fast computation, it will not make a new copy, 
but simply store the address. Since no member function in the class modifies the 
superstate data, an application using or generating such an array for other purposes 

30 may expose it to a Quantifier object 600 without risk.) Where multiple method 
options exist, as in the choice of how to perform integration or volume estimation, 
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these are specified by flags 614. Member-setting functions 640 may modify either 
the default or particular values of these flags. 

Without valid members of the above types, the Quantifier object 600 is not a 
valid object. Other data members are needed only for certain functions and may be 

5 assigned as necessary. There are optionally assigned arrays 615 for bound density 
and 616 for free density. The arrays 617 for binding rate and 618 for creation rate 
are most often used in the context of sequential data, but are available for a single 
Quantifier object 600. If segmentation software has classified the points in the grid 
according to enumerated classes si, s 2 , . . then an application may make an IxJxK 

10 array 619 of corresponding classifiers available to the Quantifier object 600, 

communicating its address with a member-setting function 630. If another function 
in the application has constructed triangulated surfaces for any of the resulting 
segments, or other surfaces of interest, the function may provide a list 620 of 
pointers to such surfaces. 

1 5 The function of the invention may also provide a pointer 621 to a transport 

model, such as those discussed in (copending U.S Serial No. 09/566,478) which 
transport model is applied to a function in a standardized format to give the 
predicted vector of flow at a point (represented by the example of D • Vc in an above 
equation, so that the Quantifier can compute flow through a surface without itself 

20 containing a fixed transport model. Similarly, solving an equation such as (16) for 
c, given the right hand side, is a service to be provided by the transport model, 
which is equipped with the appropriate explicit form for the left hand side. The 
quantifier object does not and need not "know" whether the transport consists of 
pure diffusion or includes net fluid flow, or the values of such coefficients as D. 

25 The function members 630 of a Quantifier object 600 act on the data it 

contains. The volume-in-range functions 631 computes the volume of the region in 
which the superstate density c lies between limits C 0 and Ci, within the domain of 
the data as default, or within a specified segment or surface, from the list 61 1 or 
elsewhere. Where the default lower limit 0 or default upper limit oo is used, which a 

30 density value cannot lie outside, the associated computation may be skipped. 
Optional arguments, taking default values if unspecified, are shown in square 
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brackets [ ]. The function performs the computation by one of the methods shown 
in Drawing 8, or by another method added within the architecture of the present 
invention, according to the method flag 614 currently in place. The integration 
function 632 returns an estimate of the total quantity of a superstate with density c g9 

5 within the domain of the data or within a specified segment or surface, from the list 
620 or elsewhere. The function performs the computation by one of the methods 
shown in Drawing 9, or by another method added within the architecture of the 
present invention, according to the method flag 614 currently in place. Given as 
argument the address of a triangulated surface, either previously declared segment 

1 0 surface or an arbitrary surface in the same format, the surface integration function 
633 computes the flux through the surface, as shown in Drawing 10. The difference 
function 634 computes the point-wise difference between the data array 613 of the 
Quantifier 600 and the subtrahend 635 to be subtracted from it, with optional 
arguments that can restrict the computation to a particular segment, specify a divisor 

1 5 (such as the size of a time step) to applied to the difference before storing it ? and 
give a target array (rather than an automatically created new array) or a complete 
Quantifier object 600 where the results are to be stored. A generalized version 636 
allows the application programmer to specify a list of arrays and the code address of 
an algorithm to be applied to their entries, for such purposes as constructing a 

20 Kalman filter to be applied to a sequence of data arrays. An image processing 

function 637 allows specification of the code address of a local operator such as a 
discrete Laplacian or an edge detector; a library of commonly-used operators is 
included, and may be referenced by name. Since an edge detector may return 
different values of edge strength in different directions, the library and API include 

25 a vector-array class to receive such values. 

Drawing 7 shows the Temporal Quantifier object 700 that handles a 
temporal sequence of volume data sets, whether arising by scan or by simulation. 
To do so, its data members 710 include sequences of members corresponding to 
single members of a static Quantifier object 600. It also contains corresponding 

30 timeless members: array dimensions 711, scalar type 714, method choice flags 716 
and transport model 722. Being identical to the corresponding members of a static 
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Quantifier 600, in one embodiment these timeless members are inherited from a 
common "abstract" C++ class. In a Java implementation of the present invention, 
however, the Quantifier and Temporal Quantifier objects 600 and 700 would be 
coded separately. Additional members specific to the temporal case are the number 

5 712 of instants for which datasets (in a common array size and scalar type) are to be 
analyzed, and the times 713 of those instants. The data arrays 715 are now a list, 
one for each of the times 713. If segmentation 721 is static, with the substrate 
segmented in the same way throughout an analysis, it may be represented by a 
single array: if it is dynamic (as is necessary, for example, in the study of diffiision 

1 0 through a moving lung), it must be represented by arrays for each time. 

By default, the volume, integration and efflux functions 731, 732 and 733 
respectively produce the sequences of lvalues that would result from applying 
those in the static Quantifier 600 to the N arrays in sequence: an optional argument, 
applicable also to the other function members 730 of a Temporal Quantifier object 

15 700, allows selection of a restricted range or subset of times. These functions yield 
sequences of numbers: further functions 734, 736 and 737 store results in sequences 
of target arrays, creating created new array sequences for this purpose if necessary, 
or in complete Temporal Quantifier objects 700, as specified by the way the 
function is called. The rate of change function 734 makes use of the subtraction 

20 function 634, and its generalization 736 correspondingly exploits the multi-image 

function 626 in a sequential manner. The operator 738 in the sequential version 737 
of the image processing function 637 may refer to data points by spatial offset from 
a current point within the current data array, by temporal offset in the sequence of 
arrays, or by a combination of such offsets. 

25 A group 740 of member-setting functions allow fixing or changing of the 

various data members of the Temporal Quantifier object 700. 

Drawing 8 shows the logic which implements formula (10) for the volume 
contained in a triangulated surface, on the assumption that the surface is 
topologically closed (each edge is shared by exactly two triangles). The function 

30 implementing this logic receives a triangle list as input, and initializes 800 an 

iteration procedure to the first triangle and a floating point number Volume 6 to 0.0. 
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It reads the (x,y,z) coordinates of the current triangle a, as in (9), evaluates 810 the 
expression 

((xM- x 0 (a))(y 2 (a)-y Q (a))- (x» - xMlyfc)- )))(^) + + 
adding it to Volume 6, and checks 820 whether this is the last triangle. If No, it 
5 goes 821 to the next triangle given by the iterator; and repeats step 810; if Yes, it 
returns 850 the accumulated Volume6, divided by the common factor 6.0, as the 
contained volume. 

A coarse estimate from the sampling c g of the volume for which c is between 
C x and C 2 follows the flow shown in Drawing 9. The function implementing this 

1 0 logic receives array details, the values of C\ and C 2 , and optionally the details of a 
segment S. It sets up 900 an iterator over points that may be in 5, given the 
information (bounding box, octree, etc.) in the segmentation data, and initializes an 
integer counter to 0. It tests 910 the current point for lying between d and C 2 , and 
optionally 91 1 for being a point of the segment S. (The order of these steps may be 

1 5 reversed, depending on efficiency considerations depending on the architecture of 
the central processor and on the fraction of not-in-S points that may be presented by 
the iterator.) If the point is in range and (optionally) in S, the count is increased 912 
by 1 . The iterator then tests 920 whether this is the last point: if not, it moves 921 
to the next, and repeats the tests 910 and/or 911. If there are no more points, it 

20 returns 950 the value of the count, multiplied by the volume H=h x h y h z of the unit 
volume element. A similar logic flow estimates the quantity of superstrate 
represented by a density c (x,y,z) over a segment 5, replacing the integer count 
increment 912 by addition of the floating point number c (x,y,z) or a more 
sophisticated (and numerically costly) estimate of the local contribution to the total. 

25 The preferred embodiment of the present invention uses simple addition at this 
point. 

A more accurate volume estimate follows the logic shown in Drawing 10. 
For simplicity, we show here the version that computes the volume of the region 
where c (x,y,z) is above (or below, according to a flag set as the ftmction is called) a 
30 single threshold C. For the volume between C x and C 2 one may either separately 
compute the volume below each value and subtract the first from the second, as 
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detailed in Drawing 12, or (where the time of memory access to obtain values of c 
(x,y,z) is important) duplicate the three steps 1010, 1020 and 1030 below, either 
adding their difference to the volume accumulator For maintaining two such 
accumulators and subtracting one from the other after all volume elements have 

5 been processed. The function initializes 1000 an iterator on the first 8-cornered 
volume element (note that an IxJxK grid has UK points but only (I-\)(J-l)(K-\) 
such elements), and sets the floating point number Vto 0.0. In step 1010 it 
determines which case T of the 256 cases shown (up to rotation and — where 
permitted — reversal) in Drawing 5 matches the combination of in-range and out- 

1 0 of-range corners of the current volume element. Representing the state of each by a 
0 or 1, in a determined sequence, this combination corresponds to an 8-bit number, 
which is used to index a look-up table to select the code used in step 1020 to 
compute the linear or other interpolation the edge lengths between in-range and out- 
of-range corners appropriate for this case, and the polynomial from (1 1) in which 

15 step 1030 inserts them to evaluate the contribution that it adds to the volume V. 
Step 1040 of the iteration determines whether a next element must be found 1041 
passed to step 1010; if not, step 1050 multiplies the accumulated Fby H-h x h y h z md 
returns the result. 

Where the required quantity is the volume within a segment S for which the 
20 value of c (x,y,z) is above C, the logic is slightly more complex. The function 

computing this quantity receives the details of S and the value of C, and 1 100 sets 
up an iterator over volume elements that may be in S 9 also initializing a volume 
estimate to 0.0. It tests 1 105 whether all corners of the current volume element are 
outside S 9 and if so passing control to the iterator test 1 140 of whether any more 
25 elements remain. It then tests 1 106 whether all corners are inside S; if so, control 
passes to a sequence 1 1 10, 1 120, 1 130 of steps identical to 1010, 1020 and 1030 
respectively. If not, it computes 1 107 a fractional contribution by a method such as 
one of those described above. 

As in the case of the volume computation without reference to a segment S, 
30 where the volume of points such that C\ <c(x,y,z)<C 2 is required we may use the 
iterator over volume elements once, computing separately the volumes below C\ 
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and below C 2 for each element and adding them to distinct running totals, 
subtracting one from the other at the end, or adding their difference to a single total. 
Alternatively, as in Drawing 12, we may 1200 input the data, independently set 
range conditions "<Ci" and M <C 2 " in steps 1210 and 1211 respectively, apply 1220 
and 1221 the whole process FV and FV t represented respectively in Drawing 10 or 
1 1 with these inputs, and 1250 report the difference of the two results. 

Given grid values of c at two or more successive time points, numerous 
methods exist for estimating dc/dt, of which the simplest is to subtract the most 
recent value at each point from the current value there, and divide the result by the 
time difference. This is equivalent to a straight line fit at each grid point. More 
complex methods allow fitting of more complex curves to a longer sequence of 
values, and smoothing across neighboring points. In the present invention we 
choose one of these available methods and implement it both as a callable function 
in the API and as step 1310 of the logical flow shown in Drawing 13 of a function 
which estimates the unbound portion J of a measured density c g . The function 
stores this estimate in an array with address Z), then selects 1320 a boundary B 
within which are above-noise values of c but on which c is negligible or (if no such 
B exists) throws an exception. It fixes 1330 the address A of an array a (creating the 
array if necessary, or accepting an address passed to it) in which values of<f are to 
be stored. It passes 1340 the data D, A and B to a transport model function which is 
required to fill a with estimated values of <f by inverting the transport operator, and 
upon a successful return from this function returns 1350 the address A. 

Another ftmction in the present invention that requires cooperation by the 
transport model is the computation of flux through a triangulated surface. Its input 
is a list of N triangles, typically but not necessarily those in a closed surface such as 
the boundary of a segment, with the vertices of each listed in anticlockwise order as 
seen from outside. As shown in Drawing 14, it begins 1400 by creating an array e 
of N+l scalar values, initialized to 0.0. It initializes (step 1410) to 0 the index I of 
the current triangle in the input list and the current scalar in the array e. It defines 
1420 edge vectors u\ and u 2 by subtracting the coordinates of the zeroth vertex of 
the current triangle from those of the first and second respectively, and computes 
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1430 the mid-point w of the current triangle. It then 1440 passes w to the transport 
model, receiving 1441 the flow vector F returned by the transport model for that 
point. It computes 1450 the estimate w • F of flux through triangle /, adding it to the 
final number in the array e (which does not correspond to an individual triangle) and 
storing the resulting per-area flux in e[i\ (for uses such as display). It tests 1460 
whether there is another triangle in the list, if so advancing 1461 to the next triangle 
(by incrementing I) and returning to step 1420, otherwise returning the address of e. 
Alternatively, if the calling function creates the array e and supplies its address as 
input, the function here described may return only a signal of successful completion. 

The invention is particularly useful within medical procedures for 
determination of a specific administration therapy that should be provided. A 
volume of the body of the patient is volumetrically evaluated, the available locations 
or a specific location is assumed, the dynamic response (migration, movement, 
dispersal, absorption, adsorption, etc.) of administration of a material (e.g., drug, 
marker, toxin, cell, etc.) at a specific point or at the various points is interpreted or 
estimated or quantified (hereinafter collectively referred to as "determined") on the 
basis of the volumetry analysis or evaluation, and a therapy is selected on the basis 
of results of the interpretation, estimation or quantification. That determined 
therapy may then be considered for use, and ultimately used upon approval of the 
determined therapy or procedure. 



27 



